import numpy as np
import matplotlib.pyplot as plt


def hz2mel(hz):
    return 2595. * np.log10(1 + hz / 700.)


def mel2hz(mel):
    return (10. ** (mel / 2595.) - 1.) * 700.


sr = 16000.
n_spec = 257
n_filters = 40
low_hz = 0.
high_hz = sr / 2.
low_mel = 0.
high_mel = hz2mel(high_hz)
mel_arr = np.linspace(low_mel, high_mel, n_filters + 2)
hz_arr = mel2hz(mel_arr)
bins = np.floor(hz_arr / high_hz * n_spec)

spr = 2
spc = 1
spn = 0
plt.figure(figsize=[8, 8])
spn += 1
plt.subplot(spr, spc, spn)
plt.title('bins')
y1y2 = np.tile([[100], [-100]], len(bins))
plt.plot([bins, bins], y1y2)

filters = np.zeros((n_filters, n_spec), dtype=np.float64)
for i in range(1, n_filters + 1):
    m_minus = int(bins[i - 1])
    m_center = int(bins[i])
    m_plus = int(bins[i + 1])
    for k in range(m_minus, m_center):
        filters[i - 1, k] = (k - bins[i - 1]) / (bins[i] - bins[i - 1])
    for k in range(m_center, m_plus):
        filters[i - 1, k] = (bins[i + 1] - k) / (bins[i + 1] - bins[i])
spn += 1
plt.subplot(spr, spc, spn)
plt.title('bins')
for f in filters:
    plt.plot(f)

print('Please check and close the plotting window to finish ...')
plt.show()
